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OBJECTIVE 


Applied  Research  Laboratories,  The  University  of  Texas  at  Austin  (ARL:UT),  proposes 
to  incorporate  digisonde  data  into  a  global  ionospheric  analysis  algorithm  to  provide 
more  accurate,  global  specifications  of  the  ionospheric  space  weather.  This  work  will 
improve  the  Ionospheric  Data  Assimilation  Three  Dimensional  (IDA3D)  algorithm, 
which  produces  global  electron  density  specifications.  The  IDA3D  will  be  modified  to 
incorporate  the  return  time  spectrum  (or  virtual  height  profile)  observed  by  digisondes 
distributed  around  the  world.  The  improved  specifications  can  be  made  available  to  the 
United  States  Air  Force  (USAF)  operational  planners,  allowing  them  to  respond 
appropriately  to  the  ionospheric  weather. 

WORK  PREFORMED 


The  tasks  preformed  under  this  grant  falls  into  three  major  categories:  development  of  the 
forward  model  H,  modification  of  IDA3D  to  incorporate  nonlinear  observations,  and 
validation  of  the  new  IDA3D  results.  In  addition,  it  was  necessary  to  become  familiar 
with  the  digisonde  observations  and  data  availability.  To  become  familiar  with  the 
digisonde  data,  two  days  of  data  from  the  Digital  Ionogram  DataBase 
(http://ulcar.uml.edu/DIDBase/).  the  largest  digisonde  database,  were  visually  examined. 
These  data  were  processed  by  the  Automatic  Real  Time  Ionogram  Scaler  with  True 
Height  (ARTIST)  [Reinisch  and  Huang,  1983]  program  into  electron  density  profiles, 
and  the  IDA3D  algorithm  constructed  space  weather  maps  for  the  29-30  October  2003 
magnetic  storm  with  these  profiles.  The  processing  of  the  digisonde  data  for  this  storm 
study  provided  an  extensive  introduction  into  the  data  timeliness,  quality  control,  and 
database  issues  involved  in  using  digisonde  data. 

The  first  task  was  to  develop  a  forward  model  to  predict  the  virtual  height  from  a  given 
array  of  gridded  electron  densities.  The  forward  model  was  developed  from  the  theory 
developed  in  Budden  [1961].  The  virtual  height  h’(f)  at  a  given  digisonde  emission 
frequency /is  given  by 

*'(/)=  ]>(/./,(*))*  (>) 

0 

where  ju’  is  the  real  component  of  the  group  index  of  refraction,/,  is  the  electron  plasma 
frequency,  and  zr  is  the  reflection  height  of  the  wave.  The  electron  plasma  frequency  is 
related  to  the  electron  density  by  fp=c^n  where  n  is  the  electron  density  and 
c=8. 97866275  Hz  m'1 5.  The  reflection  height  zr  is  given  by  fp{zr)  =  f  The  real  group 
index  of  refraction  is  given  by 
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where  B  is  the  strength  of  the  magnetic  field  at  a  given  geographic  latitude  X,  geographic 
longitude  <j>,  and  altitude  z,  ©  is  the  angle  between  the  magnetic  field  line  and  the 
propagation  direction  of  the  wave  (assumed  to  be  vertical),  X  =  (fp/f)2  is  the  ratio  of  the 
plasma  frequency  to  the  emission  frequency,  Y=gpB/f  is  the  ratio  of  the  electron 
gyrofrequency  to  the  emission  frequency  and  gp=2. 799249247  x  1010  C/kg,  the  subscripts 
L  and  T  denote  the  longitudinal  (cos(©))  and  transverse  (sin(©))  components  of  the 
waves,  m  is  the  mode  factor  with  m=  1  for  the  ordinary  (0)  mode  of  the  wave  and  -1  for 
the  extraordinary  (X)  mode,  //  is  the  real  phase  index  of  refraction  given  by 


V"~A  ' 
*  V  D 


(3) 


and  D  and  a  are  simplifying  algebraic  functions  given  by 


D  =  \-X-\Y2  +m4ct 


and 


+yL2(i-xy 


As  X  approaches  1,  ji'  goes  to  infinity,  and  fj.  goes  to  0.  However,  Budden  [1961] 
published  a  reliable  approximation  to  fi '  for  both  the  O-mode  and  the  X-mode  when  is 
near  the  reflection  altitude.  The  O-mode  corresponds  to  reflection  at  X=\.  The 
approximation  that  is  valid  for  1-A«1  is 


A, 


1 


approx 


\  3(1  -X) 3 

2 tan2  © 


(4) 


sin  ©Vl  -  X 

The  X-mode  corresponds  to  reflection  at  X=l±Y  with  an  approximation  for  \±Y-X«\ 
produces 


2  +  7 


* approx 


(5) 


^2(1  +  7  -  X)(\  +  7)(1  +  cos2  ©) 
where  the  negative  sign  corresponds  to  the  X=l-Y  reflection  and  the  plus  sign 
corresponds  to  th eX=l+Y  reflection.  The  forward  model  for  a  given  electron 

density  distribution  is  a  numerical  approximation  to  the  integral.  The  electron  densities 
specified  on  a  grid  of  L  horizontal  points  (irregularly  spaced  in  latitude  and  longitude) 
and  K  vertical  points  with  the  densities  placed  in  the  center  of  the  grid  cell.  The  forward 
model  is  given  as 


H(nl<k)  =  Y,nXnin)dz+  j  p\ 


approx 


dz 

WP 


dz 


(6) 


The  first  term  in  the  forward  model  is  a  loop  through  the  grid  over  those  frequencies  less 
than  the  emission  frequency.  The  electron  density  is  linearly  interpolated  between  the 
surrounding  data  points,  and  the  vertical  step  distance  dz  is  step  to  1.0  km/p’.  As  the 
forward  model  steps  in  altitude,  the  approximation  (4)  is  compared  with  the  numerical 
solution  (2).  When  the  two  values  are  within  5%,  the  integral  is  solved  analytically.  For 
the  O-mode,  the  analytic  term  reduces  to 
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approx 


dz  ,  dz 
-dz  = 
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4  tan  0 


arccos(l  -  2(1  -  X ))  + 


3^X(l-X) 
2tan2  © 


(7) 


The  altitude  derivative  is  the  inverse  of  the  slope  of  the  linear  interpolation. 


This  forward  model  is  dependent  upon  both  a  background  electron  density  and  a 
magnetic  field  model.  An  ideal  3DVAR  algorithm  does  not  utilize  any  inputs,  which  are 
not  within  the  background  model.  However,  the  magnetic  field  within  the  ionosphere  is 
only  weakly  dependent  upon  the  plasma  conditions.  During  magnetic  storms,  the  changes 
in  the  low  altitude  magnetic  field  strength  are  less  than  a  few  percent.  In  addition,  the 
magnetic  field  is  a  secondary  term  for  the  O-mode  virtual  heights.  Hence,  the  magnetic 
field  is  considered  to  be  static  in  time,  and  an  independent  magnetic  field  model  does  not 
violate  the  spirit  or  the  functionality  of  3D  VAR.  At  present,  the  International 
Geomagnetic  Reference  Field  (IGRF)  model  was  used  for  the  magnetic  field.  The 
forward  model  adequately  calculates  the  virtual  height  distribution  from  a  given  electron 
density  profile. 


The  second  task  was  to  implement  the  forward  model  into  IDA3D.  The  fundamental 


xa  =x6  +  PbHT[HPbHT  +P j\yo-H(&))  (8) 


3DVAR  equation,  which  IDA3D  solves,  is 

where  xa  is  the  analysis  vector,  xb  is  the  background  model  vector,  Pb  is  the  error 
covariance  matrix  associated  with  the  background  model,  P0  is  the  error  covariance 
matrix  associated  with  the  observations,  H  is  the  transformation  matrix,  y_0  is  the 
observation  vector,  and  H  is  the  forward  model.  For  a  linear  forward  model,  the 
transformation  matrix  is  equal  to  the  component  parts  of  the  forward  model,  H(x)  =  Hx. 
However,  the  transformation  matrix  is  different  for  a  nonlinear  forward  model.  It  is  given 
by  the  derivative  with  respect  to  x  at  the  background  model.  For  the  virtual  heights,  the 
elements  of  the  transformation  matrix  are  given  by 


H,m  = 


dXh\ 


dn . 


dfp  dn 


b,m 


(9a) 


outside  of  the  analytic  regime 


/  dz 

fl-  3  1 

2  ,  3(l  -2  A) 

8X 

2 sin©  dfp 

v.  4tan2©; 

yl2(\-X)  4  tan 2  ©-JX(\-  X) 

dnb,m 

(9b) 


and  inside  of  the  analytic  regime.  The  subscripts  l  and  m  correspond  to  the  observation 
and  grid  elements,  respectively.  It  should  be  noted  that  the  second  two  terms  in  equation 
(9)  have  contributions  to  the  grid  cells  above  {m+ 1)  and  below  (m- 1)  the  reflection 
height. 
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In  addition  to  developing  the  transformation  matrix,  the  basic  convergence  loop  of  the 
algorithm  was  modified.  Because  the  matrix  inversion  in  equation  8  is  numerically  and 
computationally  difficult  and  only  the  solution  is  important,  equation  8  is  solved  by  a 
steepest  descent  algorithm,  which  iterates  to  convergence.  After  the  first  iteration,  the 
previous  iteration’s  analysis  vector  is  used  as  the  background  vector.  For  the  nonlinear 
observations,  the  transformation  matrix  is  updated  at  every  iteration,  and  the  combined 
matrix  HPbHT+P0  is  recalculated.  This  increases  the  IDA3D  runtime,  but  not 
prohibitively.  This  new  convergence  loop  allows  IDA3D  to  treat  nonlinear  data. 


The  last  task  was  to  conduct  test  of  the  new  data  source.  The  test  period  for  this  study 
was  October  30,  2003.  While  an  extreme  storm  period,  it  has  been  previously  studied 
with  IDA3D.  The  data  is  locally  available  and  previously  quality  checked.  In  addition, 
IDA3D  maps  using  ARTIST-calculated  profiles  from  hand  scaled  ionograms  are 
available  for  comparison.  The  first  test  run  of  the  EDA3D  used  only  O-mode  autoscaled 
virtual  height  profiles  from  five  different  digisondes.  The  assumed  error  in  the  virtual 
height  was  10%.  The  background  model  used  was  the  Riley-ICED-Bent-Gallagher 
(RIBG)  model.  Figure  1  shows  a  comparison  between  the  measured  digisonde  virtual 
heights  and  the  IDA3D  and  RIBG  calculated  values  for  the  Osan  digisonde  (URSI  code 
SN437,  37.1°  lat,  127.0°  long),  while  Fig.  2  is  a  similar  comparison  between  the  IDA3D, 
RIBG,  and  ARTIST  calculated  electron  density  profiles.  Figures  3  and  4  are  similar  plots 
for  the  Dyess  digisonde  (URSI  code  DS932,  32.4°  lat.,  260.0°  long).  A  similar  run  was 
also  conducted  with  the  full  set  of  available  observations,  including  GPS  slant  total 
electron  content  (TEC)  from  ground  and  satellite-based  receivers  and  in  situ  electron 
density  measurements  from  orbiting  spacecraft.  As  these  figures  show,  IDA3D  is  able  to 
reproduce  the  measured  virtual  height  profiles  well.  The  IDA3D  electron  density  profiles 
produced  by  including  the  virtual  heights  are  also  in  good  agreement  with  the  ARTIST 
profiles.  One  interesting  feature  of  this  study  is  the  development  of  the  valley  region  in 
the  IDA3D  analysis.  These  valleys  tend  to  be  thin  layers. 


Osan 


Frequency  (MHz) 

Figure  1.  A  comparison  between 
the  measured  virtual  height  profile 
(solid  line)  with  a  virtual  height 
profile  calculated  from  the  RIBG 
model  (dashed  line)  and  the  IDA3D 
analysis  (diamond  symbols)  for  the 
Osan  (URSI  SN437)  digisonde. 


Oson 


Figure  2.  A  comparison  between 
ARTIST  produced  electron  density 
profiles  (+)  with  the  RIBG  model 
(dashed  line)  and  the  IDA3D 
analysis  (diamond)  for  the  Osan 
station. 
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Dyess 


Figure  3.  A  virtual  height  comparison 
for  the  Dyess  digisonde.  The  solid  line 
shows  the  measured  values,  the  dashed 
lines  gives  the  profile  calculated  from 
RIBG,  and  the  diamonds  show  the 
IDA3D  virtual  heights. 


Dyess 


Figure  4.  A  comparison  between 
ARTIST  produced  electron  density 
profiles  (+)  with  the  RIBG  model 
(dashed  line)  and  the  IDA3D 
analysis  (diamond)  for  the  Dyess 
station. 
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NEW  FINDINGS 

The  measured  virtual  height  traces  can  be  directly  ingested  into  a  3DVAR  algorithm.  In 
doing  so,  reasonable  agreement  is  reached  between  the  3DVAR  solution  and  the  solution 
presented  produced  from  the  state-of-the-art  inversion  model.  In  addition,  the  analysis 
profile  produces  the  valley  region  between  the  E  and  F  layers. 
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